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Abstract: We study the 9 dependence of four-dimensional SU(iV) gauge theories, 
for N > 3 and in the large- N limit. We use numerical simulations of the Wilson 
lattice formulation of gauge theories to compute the first few terms of the expansion 
of the ground-state energy F(9) around = 0, F(9) - F(0) = A 2 9 2 (1 + b 2 9 2 + ...). 
Our results support Witten's conjecture: F{9) — F(0) = A9 2 + 0(1 /N) for suffi- 
ciently small values of 9, 9 < ir. Indeed we verify that the topological susceptibility 
has a nonzero large- iV limit Xoo = 2^4 with corrections of 0(1/ N 2 ), in substan- 
tial agreement with the Witten-Veneziano formula which relates Xoo to the r/ mass. 
Furthermore, higher order terms in 9 are suppressed; in particular, the 0(9 4 ) term 
b 2 (related to the rj' — rf elastic scattering amplitude) turns out to be quite small: 
&2 = —0.023(7) for iV = 3, and its absolute value decreases with increasing N, 
consistently with the expectation b 2 = 0(1/N 2 ). 



Keywords: 



Gauge Field Theories, 1/N Expansion, Lattice Gauge Field Theories . 



Contents 

|l|. Introduction [l] 

|2|. Numerical simulations |3| 

g. Results § 



1. Introduction 

It is rather well established that four-dimensional SU(iV) gauge theories must have 
a nontrivial dependence on the angle 9 that appears in the Euclidean Lagrangian as 

C e = \F^x)F^x)-i9q(x) (1.1) 

where q(x) is the topological charge density 

q{x) = ^W^>)^)- (1-2) 

Indeed, the most plausible explanation of how the solution of the so-called U(1)a 
problem can be compatible with the 1/N expansion (performed keeping g 2 N fixed 
[0]) suggests that in the pure gauge theory the ground-state energy depends on 9 to 
leading order in 1/N [[|. Numerical nonperturbative studies, based on Monte Carlo 
simulations of the lattice formulation of SU(iV) gauge theories (see e.g. Ref. || |J 
for recent reviews), have given nonzero estimates for the topological susceptibility x 
at 9 = 0, which is the zero-momentum two-point correlation function of q(x), and 
therefore the second derivative of the ground state energy with respect to 9 computed 
at 9 = 0. 

The 9 dependence of the theory is described by the partition function 

Z{9) = J [dA] exp (- J d 4 xC e ^j = exp[-VF{6)] (1.3) 

where F{9) is the free energy (ground state energy) per unit space-time volume. 
Some information on the behavior of F{9) comes from large-iV arguments. Witten 
argued that in the large- iV limit F{9) is a multibranched function of the type 

F{9) = N*mm k H (^^) (1-4) 
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which is periodic in 9, but not smooth since at some value of 9 there is a jump 
between two different branches. More recently, the conjecture was refined || leading 
to a rather simple expression for AF(9) = F(9) — F(0) in the large- iV limit, that is 

AF(9) = A min k (9 + 2nk) 2 + 0(1/N). (1.5) 

In particular, for sufficiently small values of 9, i.e. 9 < ir, 

AF{9)=A9 2 + 0{1/N). (1.6) 

Thus possible 0(9 4 ) terms are expected to be depressed by powers of 1/N. This 
conjecture has been supported using arguments based on a duality of large- N gauge 
theories with string theory |J . It has also been discussed in a field-theoretical frame- 
work in Ref. 0. 

In this work we study the 9 dependence of SU(iV) gauge theories, using numerical 
simulations of their Wilson lattice formulation. Numerical Monte Carlo studies of 
the ^-dependence are made very difficult by the complex nature of the 9 term. In fact 
the lattice action corresponding to the Lagrangian ( |1 . 1|) cannot be directly simulated 
for 9 0, not to mention the difficulty to define a topological charge density on the 
lattice. Here we restrict ourselves to the region of small 9 values, where one may 
expand F{9) around 9 = 

AF(9) = J2^9 2t . (1.7) 

i=l 

The coefficient A 2 gives the topological susceptibility, i.e. 

A 2 = \x, X = X2 = ^(Q 2 )e=o, (1.8) 

where Q is the topological charge Q = J d 4 xq(x). Higher order coefficients of the 
expansion ( |1 . 7| ) can be related to higher moments of the probability distribution p(Q) 
of the topological charge. For instance, 

Ai = ~2l X4 ' 



x, = v 



<QVo-3((Q 2 ) e= o) 2 • (1.9) 



A convenient parametrization of the free energy is achieved by introducing the di- 
mensionless function 

f(9) = a- 2 AF(9), (1.10) 

where a is the string tension at 9 = 0. The coefficients of the ^-expansion of f(9) 
are dimensionless scaling quantities. Let us write f(9) in the form 

f(9) = 1 -C9 2 s(9) (1.11) 
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where 



^ = ^ = 4 (1-12) 
s(9) = l + b 2 9 2 + b 4 9 4 + ..., (1.13) 

& 2 = 4 1 = , (1-14) 
A 2 12*' 1 J 

The function s(#) parametrizes the corrections to a simple Gaussian behavior. Wit- 
ten's conjecture implies that C has a finite nonzero large- N limit with 0(1/N 2 ) 
corrections, and b 2 i = 0{1/N 2 ). In particular one expects that b 2 is small, and that 
it should rapidly decrease with increasing N. We recall that the nonzero value of Xoo 
is essential to provide an explanation to the U(1)a problem, and can be related to 
the rj' mass |2|, || through the relation 

f 2 TYb 

Xoo = 2A= J -^ + 0(l/N). (1.15) 

The quantity b 2 also lends itself to a physical interpretation, being related to the 
rj' — rj elastic scattering amplitude M. 

We present results for four-dimensional SU(iV) gauge theories with iV = 3,4,6. 
In particular, we determine the ratio C = xjo 2 and the coefficient b 2 of the 0(6 2 ) 
term in the function s(9). Once verified that the results for N = 3, 4, 6 are consistent 
with the expected behavior of the type A + B/N 2 , the large- iV limit is determined 
by a fit. This allows us to obtain a rather accurate estimate of C^. The estimated 
values of b 2 turn out to be very small, of the order of a few per cent, and appear to 
decrease rapidly with N, supporting Witten's conjecture. Our results provide a first 
quantitative study of the function s(9), cf. Eq. ( |1.11| ), related to the corrections to 
the 9 2 term in the free energy F(9), i.e. to a simple Gaussian behavior. 

The paper is organized as follows. In Sec. ^| we describe the Monte Carlo simu- 
lations that we have performed. In Sec. |3] we present and discuss the results. 



2. Numerical simulations 

We consider the Wilson formulation of lattice gauge theories: 

S = —N(3 £ Tr [U^x)U u (x + ^{x + u)U}(x) + h.c] , (2.1) 

where U^x) G SU(iV) are link variables. In our simulations we employed the 
Cabibbo-Marinari algorithm | | to upgrade SU(iV) matrices by updating their SU(2) 
subgroups (we selected iV(iV — l)/2 subgroups and each matrix upgrading consists 
of N(N — l)/2 SU(2) updatings). This was done by alternating microcanonical over- 
relaxation and heat-bath steps, typically in a 4:1 ratio. Tables III ^ and |3] contain 
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p 




L 


iVsw/10 6 


T Q 


10 4 X 


C 


b 2 


5.9 


0.2605(14) 


16 


8.7 


29(1) 


1.549(12) 


0.0336(8) 


-0.028(12) 






12 


5.0 




1.544(7) 


0.0335(7) 


-0.021(5) 


6.0 


0.2197(12) 


16 


6.4 


77(5) 


0.741(5) 


0.0318(7) 


-0.023(7) 






12 


4.0 


70(1) 


0.728(5) 


0.0312(7) 


-0.032(3) 


6.1 


0.1876(12) 


16 


1.6 


173(4) 


0.382(6) 


0.0308(9) 


-0.032(9) 


6.2 


0.1601(10) 


20 


6.9 


540(30) 


0.1989(23) 


0.0303(8) 


-0.016(8) 



Table 1: Results of our Monte Carlo simulations for N = 3. We tabulate the coupling f3, 
the corresponding value of the string tension (taken from Ref. Q for N = 3), the lattice 
size L, the number of sweeps of the simulation N sw , the autocorrelation time tq of the 
topological charge, the topological susceptibility x-, the scaling ratio C = xl° 2 an d the 
coefficient b 2 . The error on the ratio C is computed by considering the errors on \ and a 
as independent. 
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L 


iVsw/10 6 


T Q 


10 4 X 


C 


b 2 


0.335 


0.2959(14) 


12 


2.0 




2.274(22) 


0.0297(6) 


-0.04(3) 


0.337 


0.2699(23) 


16 


2.3 


85(4) 


1.675(16) 


0.0316(11) 


-0.01(5) 


0.338 


0.2642(7) 


12 


7.9 


105(4) 


1.416(9) 


0.0291(4) 


-0.012(6) 


0.341 


0.2368(6) 


12 


6.4 


207(6) 


0.896(10) 


0.0285(4) 


-0.013(6) 


0.344 


0.2160(8) 


16 


6.8 


410(20) 


0.608(7) 


0.0279(5) 


-0.016(12) 






12 


6.1 


420(20) 


0.595(18) 


0.0273(9) 


-0.019(7) 


0.347 


0.1981(5) 


16 


10.2 


805(50) 


0.425(10) 


0.0276(7) 


-0.011(8) 



Table 2: Results of our Monte Carlo simulations for N = 4. We tabulate 7 = (3/(2N 2 ), 
the corresponding value of the string tension (taken from Ref. [pX(| ), the lattice size L (in 
most cases we combine results obtained for L 4 and L 3 x 2L lattices), the number of sweeps 
of the simulation N sw , the autocorrelation time tq of the topological charge, the topological 
susceptibility x> the scaling ratio C = xl® 2 an d the coefficient b 2 - The error on the ratio 
C is computed by considering the errors on x an d cr as independent. 



some information on our MC runs respectively for iV = 3,4,6: The coupling val- 
ues 1 (3 or 7 = f3/(2N 2 ), the corresponding value of the string tension (taken from 
Refs. (!(], [|), lattice sizes and the number of sweeps. The values of /3, 7 were cho- 



sen to lie in the weak-coupling region, i.e., beyond the first order phase transition 
in the case N = 6, and beyond the crossover region characterized by a peak of 
the specific heat for N = 3,4; see Ref. |1| for a more detailed discussion of this 
point. Computing quantities related to topology using lattice simulation techniques 
is not a simple task. In a lattice theory the fields are defined on a discretized set, 



1 In order to compare Monte Carlo results for various values of N, it is useful to introduce the 
rescaled coupling 7 = l/^g^N) — (3/(2N 2 ), which is kept fixed when taking the large-A^ limit of 
the lattice theory. 
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7 




L 


iVs W /10 6 


T Q 


10 4 X 


C 


0.342 


0.3239(8) 


12 


0.5 


250(25) 


2.55(7) 


0.0232(7) 


0.344 


0.2973(5) 


12 


0.7 


466(50) 


1.79(6) 


0.0228(8) 


0.348 


0.2535(6) 


16 


1.2 


2675(500) 


0.91(7) 


0.0220(17) 






12 


4.1 


3080(120) 


1.01(5) 


0.0244(12) 






10 


0.6 




0.99(10) 


0.0240(24) 


0.350 


0.2380(6) 


12 


0.5 


> 3500 


0.76(9) 


0.0237(28) 



Table 3: Results of our Monte Carlo simulations for N = 6. All tabulated quantities are 
as m Table |. 



therefore the topological properties are strictly trivial. One relies on the fact that 
the physical topological properties are recovered in the continuum limit. Various 
techniques have been proposed and employed to associate a topological charge to a 
lattice configuration, see e.g. the review ||. Here we employ a cooling technique. 
The topological charge has been measured on cooled configurations, using the stan- 
dard twisted double plaquette operator. As is well known, the sum over the whole 
lattice of the twisted double plaquette, measured on cooled configurations, takes on 
values Qt — ka, where k is an integer and a < 1. We determine the typical value of 
a by minimizing the average deviation of the twisted double plaquette from integer 
multiples of a, and then assign to Q the integer closest to Qt/ot. This method elim- 
inates the need for expensive, protracted cooling; usually fewer than 20 steps suffice 
(when using N(N — l)/2 subgroups) in order to observe a substantial convergence 
of the results, which become indistinguishable from results obtained using either im- 
proved operators for Q or more cooling steps. We measured the topological charge 
Q typically every 100 sweeps. 

As already observed in Ref . [ 1(J , a severe form of critical slowing down affects the 
measurement of Q, posing a serious limitation for numerical studies of the topological 
properties in the continuum limit, especially at large values of N. In Tables [l], § and|3] 
we report the estimates of the autocorrelation time tq obtained by a blocking analysis 
of the data, where the quoted errors on the estimates of tq are just indicative, and 
attempt to quantify the uncertainty in extracting tq from the blocking procedure 
outlined in Ref. ||10|| . These results are suggestive of an interesting phenomenon: As 
shown in Fig. [I], the data turn out to be well reproduced by an exponential behavior 
of the type 

In tq « N(e N £ a + c N ) (2.2) 

where £ CT is the length scale associated with the string tension, i.e. £ CT = cr^ 1 ^ 2 , 
and ejv, cn are constants. A fit to the data for N = 3,4,6, using the function 
(gj), gives respectively e 3 = 0.40(1) and c 3 = -0.42(4), e 4 = 0.42(1) and c 4 = 
—0.43(4), e 6 = 0.48(3) and c 6 = —0.55(9), showing that lnrq is approximately 
proportional to A^ (at fixed £ CT ). A qualitative explanation of this severe form of 



- 5 - 



1 

▲ 


I i | i | ;.• | i | i | i 

N=6 ik 


1 


• 


N=4 ..• 1 




1 


M Q 
IM=o 








_.■** 








1 ■ • 


'i , I 1,1,1,1, 


i 



3 I I I L_!_ I I I I I I I I I I I I 

2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 

-1/2 

a 



Figure 1: The autocorrelation time tq versus £ CT = a 1 ^ 2 . The dotted lines show the 
linear fits. 



critical slowing down may be that topological modes give rise to sizeable free-energy 
barriers separating different regions of the configuration space. As a consequence, 
the evolution in the configuration space may present a long-time relaxation due to 
transitions between different topological charge sectors. We do not have a priori 
arguments in favor of the exponential behavior ( |2.2| ), which is just a guess arising 
naturally after observing the plot in Fig. [I]. Note also that this dramatic effect 
is not observed in plaquette-plaquette or Polyakov line correlations, suggesting an 
approximate decoupling between topological modes and nontopological ones, such as 
those determining the confining properties. 



3. Results 

In this section we present and discuss the Monte Carlo results concerning the topo- 
logical susceptibility x an d the function s(9), cf. Eq. fll.llQ . In Tables |l|, || and |3 



we report \i the ratio C = xjo 2 , and the coefficient b 2 , cf. Eq. (|1.14j) . Data for 
different lattice sizes allow us to check for finite size effects: their comparison shows 
that they are under control. The ratio C is plotted in Fig. || versus a, to evidentiate 
possible scaling corrections, which are expected to be 0(a 2 ) apart from logarithms. 
Such corrections are indeed clearly observed for N = 3 and N = 4, and one needs 
to allow for them in the fitting procedure to extrapolate to the continuum limit. For 
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Figure 2: The scaling ratio C = x/® 2 versus a. The dotted lines show the linear fits. 



this purpose we assume a linear behavior 

C = C cont + ua, (3.1) 

which takes into account of the leading 0(a 2 ) scaling correction and it is also sup- 
ported by the data. We obtain 

C cont = 0.0282(12), u = 0.077(24) for N = 3, (3.2) 
C cont = 0.0257(10), u = 0.049(15) for N = 4. (3.3) 

For iV = 6 no evidence of scaling corrections is observed. Therefore, in this case as 
final estimate we take the weighted average of the data for the largest values of [3, 7, 
i.e. for 7 = 0.348 and lattice size L = 12, 16 and for 7 = 0.350, obtaining 

C cont = 0.0236(10) for N = 6. (3.4) 

Results for smaller 7 are used to control scaling corrections: they are perfectly con- 
sistent, indicating that scaling corrections are at most of the size of the error 2 . 

These results of C cont are plotted in Fig. || versus l/N 2 , which is the expected 
order of the asymptotic corrections in the large- iV limit. They are clearly consistent 

2 We mention that a fit of the N — 6 data to a constant gives C cont = 0.0232(5) with \ 2 /d.o.f. as 
0.5, while a linear fit to a + ba gives C cont = 0.0240(25) and u = -0.01(3). 



-7- 



with a behavior of the type 



Cco„t(A0 = C oo + ^. (3.5) 



Assuming such a behavior, a fit to the results for N = 3, 4, 6 gives 

Coo = 0.0221(14), B = 0.055(18), (3.6) 

which provides a rather accurate estimate of the large- N limit of C; this limit is 
clearly seen to be nonzero. Note that, extrapolating using only the N = 4, 6 results, 
we obtain a perfectly consistent result, i.e. C CO nt = 0.0219(20) and B = 0.06(4). 

In order to compare with the numbers appearing in the literature, we trans- 
late the above results in physical units assuming the rather standard value \fo = 
440 MeV, for N > 3. We obtain X o^ = 170(3) MeV, and x 1/4 = 180(2) MeV for 
N = 3. These numbers have to be compared with the r.h.s. of Eq. ( |1.15| ), which, 
using the actual values of /„., m v i, and Nf = 3, gives (191 MeV) 4 ; using the formula 
refined by Veneziano || for which m 2 , — > m 2 , + m 2 — 2m 2 K , one obtains (180 MeV) . 
Our results for N = 3 are in substantial agreement with the estimates obtained so 



far using various methods, see, e.g., Refs. [12, 13, 14, 15, 16, 17|. Ref. |9j presents an 



estimate of obtained by extrapolating results for N = 2, 3, 4, 5 using the formula 
(PI); they obtain = 0.0200(43) from C cont (N = 2) = 0.0545(25), C COIlt (N = 
3) = 0.0355(33), C cont (N = 4) = 0.0224(39) and C cont (N = 5) = 0.0224(49). Our 
work improves these results essentially for two reasons: we performed simulations 
with much higher statistics, and we do not rely on N = 2 data for our extrapolation, 
which may be too small a value of N. The large- N extrapolation of Ref. || is in 
agreement with ours, although we note some differences in the finite N results. The 
extrapolation of Ref. |J strongly depends on the N = 2 result, while their N = 2 
result is not fitted by our extrapolation (|3.6Q , which would give approximately 0.036 
for N = 2. Thus, at this stage the agreement of the large- N extrapolations appears 
to be rather accidental. The result for iV = 2 of Ref. M would suggest that N = 2 



is outside the region where the behavior ( |3.5|) sets in 



Let us now turn to the coefficient 62 in the expansion ( |1.13j ) of the function s(8), 
which determines the correction of F(9) from the simple 8 2 behavior. Note that a 
reasonably accurate computation of which is necessary to determine b 2 , requires 
very high statistics due to a large cancellation between the two terms in its definition 
(|1.9D . Sufficiently accurate results have been obtained only in the cases iV = 3, 4. 
The results are listed in Tables || and ^| for iV = 3, 4 respectively. For N = Q they were 
always consistent with zero but with a large error, essentially due to the much smaller 
statistics available, taking also into account the critical slowing down phenomenon 
discussed before. We report only the most precise result obtained by a high-statistics 
simulation: 62 = —0.01(2) for 7 = 0.348 and L = 12. The results for N = 3 and 
N = A are plotted in Fig. f|. They show scaling within the errors. Assuming that 
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0.04 




Figure 3: The scaling ratio C = x/ ' 2 versus 1/N 2 . The dotted line shows the linear fit 
to A + B/N 2 . 



the scaling corrections are negligible with respect to the typical errors of the data, 
as final estimate we consider the weighted average of the data, and as error we quote 
their typical error (in this computation we discard the data for f3 = 6.0 and L — 12 
because of possible finite size effects). We obtain 3 

b 2 = -0.023(7), for N = 3, (3.7) 
b 2 = -0.013(7), for N = 4. (3.8) 

For N = 6 we report the estimate b 2 = —0.01(2) obtained at 7 = 0.348. These results 
show that b 2 is very small and are consistent with the hypothesis b 2 = 0(1/N 2 ). One 
may also assume b 2 ~ b 2t2 /N 2 and fit the above results, obtaining 6 2 ,2 = —0.21(5). 

As already mentioned, we used a cooling technique to determine the topological 
charge of a lattice configuration. Several numerical studies have shown that cooling 
methods provide quite reliable results. Nevertheless, some of the arguments on which 
cooling relies are not rigorous. Therefore it would be worthwhile to confirm our 
results using alternative methods. A promising and more rigorous approach is the 

3 Fitting the data to a constant, one obtains b 2 = -0.0226(34) for TV = 3 and b 2 = -0.0125(36) for 
N = 4, with acceptable x 2 - These errors are substantially smaller than those shown in Eqs. ( |3.7| ) 
and ([T^). We believe that the final errors reported in the text are more reliable, since scaling 
corrections appear to be negligible only within the typical errors of the data. Nevertheless, we note 
that the 1/7V 2 behavior is verified even within the smaller errors of the fits. 
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Figure 4: The scaling ratio 62 versus a. 



one based on the index theorem for Ginsparg- Wilson fermions [|18|, using, 
e.g., the overlap lattice fermion formulation 

In conclusion, we have computed the 9 dependence of the ground-state energy 
of four-dimensional SU(iV) gauge theories around 9 = 0, by evaluating numerically 
the coefficients of its expansion to 0(9 4 ). We have verified that the coefficient of 
the 9 2 term, which is the topological susceptibility, has a nonzero large-iV limit, in 
agreement with the picture proposed by Witten [[J to provide an explanation to 
the U(l)yi problem. Corrections to the leading 9 2 term turn out to be very small. 
Indeed, the coefficient 62 of the 9 4 term is a few per cent for N = 3, and appears 
to decrease with increasing N, consistently with Witten's conjecture || implying 
62 = 0(1/N 2 ). Clearly, higher-order terms in the expansion in powers of 6 2 might 
not be suppressed by powers of N and therefore might spoil the simple 9 2 behaviour 
predicted by Witten for large N. The study of these terms requires the measurement 
of higher moments of the topological charge probability distribution and are beyond 
the scope of this work. Thus, for N > 3 the simple Gaussian form 



f{9) » -ce 2 



(3.9) 



is expected to provide a good approximation of the dependence on 9 of the scaling 
free energy, up to 0(9 6 ) terms, and therefore for a relatively large range of values of 
9. 
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